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Abstract 

The role of the exchange-correlation potential and the exchange-correlation kernel in 
the calculation of excitation energies from time-dependent density functional theory is 
studied. Excitation energies of the helium and beryllium atoms are calculated, both from 
the exact Kohn-Sham ground-state potential, and from two orbital-dependent approxima- 
tions. These are exact exchange and self-interaction corrected local density approximation 
(SIC-LDA), both calculated using Krieger-Li-Iafrate approximation. For the exchange- 
correlation kernels, three adiabatic approximations were tested: the local density ap- 
proximation, exact exchange, and SIC-LDA. The choice of the ground-state exchange 
correlation potential has the largest impact on the absolute position of most excitation 
energies. In particular, orbital-dependent approximate potentials result in a uniform shift 
of the transition energies to the Rydberg states. 



1 Introduction 

The Hohenberg-Kohn theorem |T[] of ground-state density functional theory (DFT) guaran- 
tees that every observable of a stationary physical system can be expressed in terms of its 
ground-state density. In principle, this is also true for the set of excited-state energies, and 
several extensions of ground-state DFT have been proposed [Q- [[L3|| . Accurate calculations 
of excitation energies, however, remain a difficult subject. Recently, some of us proposed 
a different approach to the calculation of excitation energies [|L4[], within the framework 
of time-dependent DFT (TDDFT) |l5j. The central idea is to use the fact that the linear 
density response has poles at the physical excitation energies and can be calculated from the 
response function of a noninteracting Kohn-Sham (KS) system and a frequency-dependent 
Kohn-Sham (KS) kernel. In this way, we obtain the shifts of the KS orbital differences 
(which are the poles of the KS response function) towards the true excitation energies. 
Recent applications [jl6|]- [^5j suggest that this method may become a standard tool in 
quantum chemistry. 



The success of any density functional method, however, depends on the quality of 
the functional employed. In this article, we investigate the relative importance of the 
approximations inherent in the TDDFT formalism for the calculation of discrete excitation 
energies of finite systems. This mainly concerns the role of the ground-state XC potential, 
v xc (r) compared to the dynamical XC kernel, / xc (r, r'; u>). For the helium and beryllium 
atoms, we compare the results obtained from using the exact XC potentials and two orbital 
dependent potentials, one based on the exact exchange expression and the other on the 
self-interaction corrected local density approximation [^6|], evaluated with the method of 
Krieger, Li and lafrate (KLI) []27|]- []3l|] in combination with three distinct approximations 
for the XC kernels, which are given in Sect. |2~2. 



2 Formalism 

2.1 Kohn-Sham equations for the frequency-dependent linear density re- 
sponse 

The frequency dependent linear density response n lcr (r, to) of electrons with spin cr, reacting 
to a perturbation v\ a < of frequency cu can be written in terms of the interacting density- 
density response function x<j<j' by [01 



„( r ,.) = £/dV w(r , r VK„,(r>). 



In the spin-dependent version [^3] of time-dependent DFT |15[| , the density response n lcr 
can be expressed in terms of the response function Xsao' of the non-interacting Kohn-Sham 
(KS) system [0, 0: 

n la (r,uj) =J2 Jd 3 r' Xsaa>(r,r';uj)v Stl(7 ,(r',Lj) . (2) 

cr' 

The KS response function 

Xsaa'{r, r;u)= 5 aa > }_^{f ka - f ja ) J - — — — (3) 

is readily expressed in terms of the unperturbed static Kohn-Sham orbitals Lp ka (with occu- 
pation numbers fj a ). Relation (^) contains the linearized KS potential 

v S!la (v,cu)=v la (v,uj) + Y: J dV + E fd 3 r'U caa ,(r,r';uj)n la ,(r',u;). (4) 

cr' I I cr' 

in which the spin-dependent exchange-correlation (XC) kernel / xc is defined as the Fourier 
transform of 

/xc ^K,^](r,t,r',tO:=^- [ ^'^ (M) 



(5) 

0T> n 0| 

Given an approximation to f xc , Eqs. (|2|) and (^) can be solved self-consistently for every 
frequency to. 
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2.2 Approximations for the exchange-correlation kernel 

For spin-unpolarized ground states, there are only two independent combinations of the spin 
components of the XC kernel, since / X cfT = /xc|j, and /xctJ, = /xcj.T : 

11 11 

/xc = ^ / ] /xc aa' = ~ (/xc TT + /xc Tl ) ^xc = ~ 2^ ° a ' f™ trer' = ^(fxc^ ~ /xc ti) ) (6) 



aa 



(contrary to common usage, we have not separated the Bohr magneton in the definition of 
G xc ). Note that f x = G x , as exchange contains only parallel spin contributions. 

The simplest possible approximation is the adiabatic local (spin-)density approximation 
(ALDA). For spin-unpolarized ground states, this leads to 



j2 horn 

/x A c LDA N(r,r') = ^(r-r^ - 



dp 2 p=n ( r)iC=0 

(7) 

where e^° m is the exchange-correlation energy per unit volume of the homogeneous elec- 
tron gas, C, is the relative spin polarization, (n-f — n{)/n, and the spin-stiffness a xc = 

&(^(p>0/p)| cmj . 

Approximate XC functionals derived from the homogeneous electron gas suffer from 
several shortcomings, such as spurious self-interaction contributions. These are very sig- 
nificant for calculations of orbital eigenvalues, as they affect the asymptotic decay of the 
ground-state potential. For example, the XC potential in the local density approximation 
decays exponentially, so rapidly that only one virtual state is bound. An alternative ap- 
proach towards the construction of improved functionals is to use perturbation theory in the 
electron-electron coupling constant^6|. This leads to orbital-dependent functionals, which 



can be solved self-consistently using the optimized effective potential method (OEP) [|37| - 
|39[| . In the time-dependent case, this method takes as a starting point a given (approx- 
imate) expression for the quantum mechanical action integral as a functional of a set of 
orbitalsffRJ. Variation with respect to a local effective potential then leads to an integral 
equation for the exchange-correlation potential. Given an exchange-correlation potential of 
that kind, the corresponding exchange-correlation kernel can be constructed in the same 



spirit [jTJ], |4l|l . The essential steps are formally identical to the OEP construction of the 



exchange-correlation potential for the ground-statep^[| 

In the time-dependent X-only approximation, A xc is replaced byfj 

A x _ only = -(1/2) £ £ [ h dt f d 3 r f d 3 r' ^(r^Mr^Mrt^rt) / |r - r'| . (8) 

a ij J -°° J J 



The orbital-dependent exchange kernel in the time-dependent KLI approximation is JR], |4l 



rTDOEP / A r 1 I Y,k fka <fka{*)V*kA T ')? , (V . 



1 In general, a Keldysh contour integral in complex time is needed ^] to avoid causality difficult iesjiit 
except when the action is local in the orbitals in time, as is the case with all approximations tested here. 
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In general, the exact x-only kernel carries a frequency-dependence. This is not accounted 
for in the present approximation However, for one- and spin-unpolarized two-electron 
systems, Eqs. ([5]) is the exact solution of the respective integral equations in the limit of a 
time-dependent X-only theory. This yields 

h=Gx= _ 2jJM^ ( i for2ete y (10) 

n{r)\r — r \n(r) \ 2|r — r | J 

Inherent to any X-only theory, the resulting kernels are lacking off-diagonal elements in 
spin space. To improve upon the X-only treatment, we use the self-interaction corrected 
(SIC) LDA @ for A xc : 

= pdt (iig> t (f),n;(f)] -E^ L c DA K(t),0] - ~£ />r /d 

00 \ ia ia 

(11) 

which is an orbital-dependent functional as well due to the explicit dependence on the orbital 
densities 

n ia (r,t) = |<Mr,t)| 2 (z = 1, 2, . . . , N/2) . (12) 

] An improvement over both ALDA and exact exchange might be provided by correcting 
ALDA for self-interaction error[^]. Within the adiabatic SIC-LDA, the exchange-correlation 
kernel reads 

/™r- SIC (r,r^) = /^(r,r',-)- 

*" -7T EJWI^MI' ( <"*fW ■ 0) + ' ) .(13) 



1 f^ r ' n ^( r ^) n U r '^y 



n OCT (r)noa(r') V Sn ka (r') 

This expression reduces to the exact result (Eq. (§)) for one electron. For more than one 
electron, spurious self-interaction parallel-spin contributions in ALDA are corrected, for both 
exchange and correlation. The correction has no affect on anti-parallel spin contributions, 
leaving simply the ALDA result. We find simply 

/£ C = / x A c LDA + A/ X S C IC Gf c c = G^ DA + A /x s c IC (14) 

where 

A s.c = _2 E. f ftgftn, 0) 1 | 

' ' " ' 1 v ' dn k (r) r - r' I 



3 Calculation of excitation energies 

The linear density response has poles at the exact excitation energies of the interacting 
system (see, e.g., [j32[). The key idea is to start from a particular KS orbital energy 
difference €j a — (at which the Kohn-Sham response function @ has a pole) and to use 
the formally exact representation (|]) of the linear density response to calculate the shifts 
of the Kohn-Sham excitation energies towards the true excitation energies Vt. To extract 
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these shifts from the density response, we cast Eq. (|^) together with (f|) into the form of 
an integral equation for n\ a : 



fd 3 y' S ai// 8(r-y')-J2 fd 3 yXsau(r,y;uj)l - + f xcu „,(y,y';u 



E J d ^yXsauir,y]uj)v lu {y,u) . 



nv(y',w) 
(16) 



In general, the true excitation energies VL are not identical with the Kohn-Sham excitation 
energies e J(T — efe CT , and the right-hand side of Eq. flTop remains finite for cu — > f2. The 
exact spin-density response n lcr , on the other hand, exhibits poles at the true excitation 
energies f2. Hence, the integral operator acting on n la on the left-hand side of Eq. fll6"|) 
cannot be invertible for to — > tt. This means that the integral operator acting on the spin- 
density vector in Eq. flip]) is non-invertible (i.e., has vanishing eigenvalues) at the physical 
excitation energies. Rigorously, the true excitation energies are those frequencies where 
the eigenvalues X(u>) of 



E / d *y' E / d 



1 



Viy - y 

A(^)7<x(r, w) 



7i + /xc W ' (y , y'; w) ) j v > (y', w) 



(17) 



satisfy 

\{Q) = 1. (18) 

For notational brevity, we use double indices q = (j, A;) to characterize an excitation 
energy; u; gcr = ej CT — denotes the excitation energy of the single-particle transition 
(ko — > jo). Consequently, we set a qa := f^ a — fj a and 



$ gCT (r) := v?L( r V>(r) 



as well as 



^(w) ;=E / rf V / My)* fr^r + /xc^(y,y'; 
i/' ^ ^ v I y y I 



(19) 



Lu)) lu ,(y',Lu). (20) 



Without any approximation, equation (|17|) can be cast [JT6]] into matrix form 



77^ - WgV + ^ 

with the matrix elements 



(21) 



M gaqlal (u) = a gV \d s r \ d A r' <&* qa (r 



J d 3 r J d 3 r' $*» (^-^ + / xc ^(r, r'; $ ? , CT ,( r ') . 



At the frequencies cu = f2, Eq. ( PT| ) can be written as 

E (Mqaq'a'(^) + ^quq'a'^qa) Pq>a>(ty = Vt(5 qa {Vt) 



(22) 



(23) 
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where we have defined 

MO) :=£ 9<J (0)/(0-av)- (24) 

The solutions f2 of the nonlinear matrix-equation (^3]) are the physical excitation energies. 
The inevitable truncation of the infinite-dimensional matrix in Eq. amounts to the 
approximation of x^ ' by a finite sum 

^ {t ,<, a) „x£c£&W. (25) 

This truncation explicitly takes into account numerous poles of the noninteracting response 
function. In any adiabatic approximation to the XC kernel, the matrix elements M qaq i a i 
are real and frequency independent. In this case the excitation energies are simply 
the eigenvalues of the (Q x Q) matrix M q(Tq / a r(Q = 0) + 5 q(T ^ a iuj qa . For bound states 
of finite systems we encounter well-separated poles in the linear density response. In our 
calculations, we include many such poles, but only those of bound states, ignoring continuum 
contributions. The nature and size of the error this introduces has been studied by van 



Gisbergen et al.[£T], and does not affect the qualitative conclusions found in this work. 

A simple and extremely instructive case is when we expand about a single KS-orbital 
energy difference uj pt 0, |Tj|. The physical excitation energies Vt are then given by the 
solution of 

A(fl) = ^ + % T ) + ... = 1. (26) 



For non-degenerate single-particle poles u pT , the coefficients in Eq. ( |26|) are given by 

A(u pT ) = Mp TpT (u pT ) (27) 



and 

t-> i \ dMp T p T 
Blu, ■ 



_| 1 \p M Vr gV (Up T )M q/(T/ pT (uJ pT ) 

Mp TpT (ujp T ) q ,^ pr U pT - Ug'v' + irj 

If the pole uo pT is p-fold degenerate, uo PlT1 = uj P2T2 = . . . = uj PpTp = lu , the lowest-order 
coefficient A in Eq. is determined by a p-dimensional matrix equation 

E^«l = W^. »=l...p, (29) 
k=i 

with p different solutions A 1 . . . A p . For excitation energies f2 close to cu , the lowest- 
order term of the above Laurent expansion will dominate the series. In this single-pole 
approximation (SPA), the excitation energies f2 satisfy Eq. reduces to 

A»(fl) « M*> = ! . (30) 
\l — Uq 



The condition fll8| ) and its complex conjugate, A*(0) = 1, finally lead to a compact expres- 



sion for the excitation energies. 

tt n ta uj + 3?A„(cj ) • (31) 
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For closed-shell systems, every Kohn-Sham orbital eigenvalue is degenerate with re- 
spect to spin, i.e. the spin multiplet structure is absent in the bare Kohn-Sham eigenvalue 
spectrum. Within the SPA, the dominant terms in the corrections to the Kohn-Sham eigen- 
values towards the true multiplet energies naturally emerge from the solution of the (2 x 2) 
eigenvalue problem 

M papa'M^ pa '(uJ ) = A£ pa (uj ) . (32) 

<t'=T4 

Then, the resulting excitation energies are: 

^1,2 = + 3ft {M pM ± M pM } . (33) 
Using the explicit form of the matrix elements (|2"2"| ) one findsQ 

fti = oo + 2^J d 3 r J d 3 r' %(r) + r '' ^ r ') ( 34 ) 

ft 2 = u + 2^J d 3 r J d 3 r' %(r)G xc (r, r'; u )Q p (r f ) . (35) 

The kernel G xc embraces the exchange and correlation effects in the Kohn-Sham equation 
for the linear response of the frequency-dependent magnetization density m(r,u) [^]. For 
unpolarized systems, the weight of the pole in the spin-summed susceptibility (both for the 
Kohn-Sham and the physical systems) at ft 2 is exactly zero, indicating that these are the 
optically forbidden transitions to triplet states. The singlet excitation energies are at fti. In 
this way, the SPA already gives rise to a spin-multiplet structure in the excitation spectrum. 
We use SPA to understand the results of different approximations, since it simply relates 
the calculated shifts from KS eigenvalues to matrix elements of the XC kernel. 

At this point we stress that the TDDFT formalism for the calculation of excitation 
energies involves three different types of approximations: 

1. In the calculation of the Kohn-Sham orbitals tpk(r) and their eigenvalues e^, one 
employs some approximation of the static XC potential v xc . 

2. Given the stationary Kohn-Sham orbitals and the ground state density, the functional 
form of the XC kernel f xcaa > needs to be approximated in order to calculate the matrix 
elements defined in Eq. 



3. Once the matrix elements are obtained, the infinite-dimensional eigenvalue problem 
( PI) (or, equivalently, (|23])) must be truncated in one way or another. 

In the following, we are going to investigate the relative importance of the approximations 
(1.) and (2.). Furthermore, truncation effects will be estimated by comparing the results 
obtained in SPA (|34| , [35| ) with the solution of the "full" problem (|3|) which is based on 



using up to 38 bound virtual orbitals. 

2 Since we are dealing with spin saturated systems, we have dropped the spin-index of $ p(7 
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Table 1: Singlet excitation energies of neutral helium, calculated from the exact XC potential 



by using approximate XC kernels (in atomic units) 









AT PiA 


(XCj 


TDOEP (x-only) 


TDOEP (SIC) 






3 




<3PA 
o r i\ 


inn a 


CO A 


lull 


CD A 


lull 


CAdtl) 


Is - 


-> 2s 


n *7/i fin 
U. / 4bU 


U.mIo 


U. lOio 


0.7687 


0.7659 


0.7674 


0.7649 


U. (D (O 


Is - 


-> 3s 


0.8392 


0.8458 


0.8461 


0.8448 


0.8450 


0.8445 


0.8448 


0.8425 


Is - 


-» 4s 


0.8688 


0.8714 


0.8719 


0.8710 


0.8713 


0.8709 


0.8712 


0.8701 


Is - 


-> 5s 


0.8819 


0.8832 


0.8835 


0.8830 


0.8832 


0.8829 


0.8832 


0.8825 


Is - 


-> 6s 


0.8888 


0.8895 


0.8898 


0.8894 


0.8896 


0.8894 


0.8895 


0.8892 


Is - 


-»• 2p 


0.7772 


0.7764 


0.7764 


0.7850 


0.7844 


0.7836 


0.7833 


0.7799 


Is - 


-» 3p 


0.8476 


0.8483 


0.8483 


0.8500 


0.8501 


0.8497 


0.8498 


0.8486 


Is - 


-» 4p 


0.8722 


0.8726 


0.8726 


0.8732 


0.8733 


0.8731 


0.8732 


0.8727 


Is - 


■» 5p 


0.8836 


0.8838 


0.8838 


0.8841 


0.8842 


0.8841 


0.8841 


0.8838 


Is - 


-> 6p 


0.8898 


0.8899 


0.8899 


0.8901 


0.8901 


0.8900 


0.8901 


0.8899 


Mean abs. dev. 


0.0022 


0.0023 


0.0021 


0.0022 


0.0020 


0.0019 


0.0017 




Mean rel. dev. d 


0.28% 


0.30% 


0.26% 


0.28% 


0.25% 


0.24% 


0.21% 





a Using the lowest 34 unoccupied orbitals of s and p symmetry, respectively. 

b Nonrelativistic variational calculation [38]. 

c Mean value of the absolute deviations from the exact values. 



4 Results for the Helium Atom 



In this section we report numerical results for excitation energies of the He atom. The 
stationary Kohn-Sham equations were solved numerically on a radial grid (i.e. without basis 
set expansion) using a large number of semi-logarithmically distributed grid points [0 up 
to a maximum radius of several hundred atomic units in order to achieve high accuracy the 
Rydberg states (n > 10) as well. 



4.1 Exact Kohn-Sham potential 

To eliminate the errors (1.) associated with the approximation for the ground-state KS 
potential, we employ the exact XC potential of the He atom to generate the stationary Kohn- 
Sham orbitals y?fc(r) and their eigenvalues e^. This isolates the effects which exclusively arise 
due to the approximations (2.) and (3.). The potential data provided by Umrigar and Gonze 



45fl were interpolated nonlinearly for r < 10 atomic units. Around r = 10 atomic units, 
the XC potential is almost identical to — 1/r. This behavior was used as an extrapolation 
of the exact exchange-correlation potential to larger distances. 

Tables |l]and ^| show the excitation energies of neutral helium calculated with the exact 
exchange-correlation potential. The results are compared with a highly accurate nonrel- 
ativistic variational calculation |^] of the eigenstates of Helium. It is a remarkable fact 
that the Kohn-Sham excitation energies Ujk = ej — are already very close to the ex- 
act spectrum, and, at the same time, are always in between the singlet and the triplet 
energiespTj |48| . 



Based on these eigenvalue differences, we have calculated the shifts towards the true 
excitation energies using several approximations for the exchange-correlation kernels f xc : 
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Table 2: Triplet excitation energies of neutral helium, calculated from the exact XC potential 
by using approximate XC kernels (in atomic units) 









AT PiA 


(xc) 


TDOEP (x-only) 


TDOEP (SIC) 






3 




<3PA 
o r i\ 


fnll a 

i un 




mil 




lull 




Is - 


-> 2s 


U. <4oU 


U.7oo7 


0. Mol 


0.7232 


0.7207 


0.7313 


0.7300 


0. 1 285 


Is - 


-> 3s 


0.8392 


0.8366 


r\ o o n o 

0.8368 


0.8337 


0.8343 


0.8353 


0.8356 


0.8350 


Is - 


-» 4s 


0.8688 


0.8678 


0.8679 


0.8667 


0.8671 


0.8673 


0.8675 


0.8672 


Is - 


-> 5s 


0.8819 


0.8814 


0.8815 


0.8808 


0.8811 


0.8811 


0.8813 


0.8811 


Is - 


-> 6s 


0.8888 


0.8885 


0.8885 


0.8882 


0.8883 


0.8884 


0.8884 


0.8883 


Is - 


-»• 2p 


0.7772 


0.7702 


0.7698 


0.7693 


0.7688 


0.7774 


0.7774 


0.7706 


Is - 


-» 2>p 


0.8476 


0.8456 


0.8457 


0.8453 


0.8453 


0.8471 


0.8471 


0.8456 


Is - 


-» 4p 


0.8722 


0.8714 


0.8715 


0.8712 


0.8713 


0.8720 


0.8720 


0.8714 


Is - 


■» 5p 


0.8836 


0.8832 


0.8832 


0.8831 


0.8831 


0.8834 


0.8835 


0.8832 


Is - 


-> 6p 


0.8898 


0.8895 


0.8895 


0.8895 


0.8895 


0.8897 


0.8897 


0.8895 


Mean abs. dev. 


0.0035 


0.0010 


0.0011 


0.0009 


0.0011 


0.0013 


0.0012 




Mean rel. dev. 


0.45% 


0.14% 


0.14% 


0.12% 


0.15% 


0.16% 


0.15% 





a Using the lowest 34 unoccupied orbitals of s and p symmetry, respectively. 

b Nonrelativistic variational calculation [38]. 

c Mean value of the absolute deviations from the exact values. 



• The adiabatic local density approximation (ALDA), with the inclusion of correlation 
contributions in the parametrization of Vosko, Wilk and Nusair [|49[| . 

• The approximate X-only time-dependent OEP (TDOEP) kernel of Eq. (|9]), which is 
based on the time-dependent Fock expression, and 

• The approximate TDOEP-SIC kernel from Eq. ( |T3D with the parametrization of Ref. 
[fPJ for the correlation contributions. 

The columns denoted by "full" show the corresponding excitation energies which are 
obtained as eigenvalues obtained from the (truncated) matrix equation (^). To investigate 
the effects of the truncation of the matrix equation (|23| ) we compare the difference between 
the single-pole approximation (SPA) and the fully coupled results. The matrix equation 
(^3|) was solved using iV = 34 unoccupied Kohn-Sham orbitals of s or p symmetry. For 
each symmetry class the resulting dimension of the (fully coupled but truncated) matrix 
in Eq. ( p3|) is (4iV x 4iV) (due to the spin-degeneracy of the KS orbitals of Helium and 
the fact that the frequency-dependent Kohn-Sham response function is symmetric in the 
complex plane with respect to the imaginary axis). Thereby, convergence of the results to 
within 10 -6 atomic units was reached within the space of bound states. When comparing 
the results from the SPA with the results from the fully coupled matrix, we observe only 
a small change in the resulting excitation energies (from a few hundredth of a percent to 
at most one half percent), independent of the functional form of the exchange-correlation 
kernel. Thus we conclude that in helium the single-pole approximation gives the dominant 
correction to the Kohn-Sham excitation spectrum. Hence, starting from the Kohn-Sham 
eigenvalue differences as zeroth order approximation to the excitation energies, the SPA can 
be used for the assignment of the excitation energies which are obtained as eigenvalues from 
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Eq. (^). Recent studies using basis set expansions [^T] indicate that further improvement 
of the fully coupled results can be expected from the inclusion of continuum states. The 
general trends of the results however, are not affected. 

In figure [I] we have plotted some typical excitation energies taken from the column 
headed "full" of table [I] and ^], We can understand the trends in this figure by analyzing 



Continuum 



0.90 



5s 
4s 



3p 

3s 



0.85 .-a 

a 

o 



0.80 



2p 



2s 



Singlet 



0.75 



Triplet 



KS 



ALDA 



TDOEP 

x-only 



TDOEP 
SIC 



Exact 



0.70 



Figure 1: Typical excitation energies of He, including the orbital eigenvalues of the exact Kohn-Sham 
potential (KS) and the corrections from time-dependent density functional theory calculated within 
the adiabatic local density approximation (ALDA), with orbital dependent functionals in the X-only 
limit (TDOEP X-only) and the self-interaction corrected version of the ALDA (TDOEP-SIC). 

the results in terms of the single-pole approximation. For the single-particle excitations 
in helium, the single-pole approximation leads to two-dimensional matrix equations for the 
excitation energies (c.f. Eqs. fl52|) - (|33])). In the following, the notation 

(6) := J d 3 r J dV$;(r)d(r,r / )$ p (r / ) (36) 
will be used for the matrix elements of the two particle operators O involved in the calcu- 
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Table 3: Singlet-triplet separations in helium obtained from Eq. (p3j), using the lowest 34 
unoccupied orbitals of s and p symmetry of the exact XC potential and employing various 
approximate XC kernels. All values are in mHartrees. 
ALDA TDOEP 



btate 


V 1 

X-only 


xc 


X-only 


sic 


exact 


2S 


42.2 


32.7 


45.2 


34.9 


29.3 


35 


11.1 


9.4 


10.8 


9.2 


7.4 


AS 


4.7 


4.0 


4.3 


3.7 


2.9 


5S 


2.4 


2.1 


2.2 


1.9 


1.4 


6S 


1.4 


1.3 


1.2 


1.1 


0.8 


2P 


16.7 


6.6 


15.6 


5.8 


9.3 


3P 


4.5 


2.6 


4.7 


2.7 


2.9 


4P 


1.8 


1.1 


2.0 


1.2 


1.3 


5P 


0.9 


0.6 


1.0 


0.7 


0.6 


6P 


0.5 


0.3 


0.6 


0.4 


0.4 


dev. c 


3.0 


1.1 


3.1 


1.3 




rel. dev. 


55% 


27% 


56% 


21% 





including correlation contributions in the form of Vosko, Wilk and Nusair |4 
b Taken from Ref. @ 
€ Mean absolute deviation from the exact values 

lation. Then, in the SPA, 

nfz lct = u p + 2(W) + 2(f xc ), nfv let = u p + 2(G xc ), AQ p = 2((W) + (f c )-(G c )), (37) 

where Ail p is the singlet-triplet splitting. Within the various approximations to the kernel, 
these levels become 

Qfz = u p +(W), {lf* lct = uj p -(W), (X-only) (38) 

= ^ + 2(W0 + 2(/ x A c LDA >, = ^ + 2(G ALDA >, ALDA (39) 

= ^+(iy) + 2(/ c ALDA )-(/ c orb ), =uj p -(W)+2(Gt LDA )-(f™ h ). SKplO) 

We begin our analysis with the splitting. In the simplest case, the TDOEP X-only 
kernel, we see that the singlet transitions are always overestimated, while the triplets are 
always underestimated. Since our TDOEP treatment is exact for exchange in this case, this 
underscores the importance of correlation. In particular, since (/ x ) = (G x ) = — (W) /2, 
the splitting is just 2(W). This matrix element is always positive, correctly putting the 
singlet above the triplet, but the splitting is typically far too big. We demonstrate the effect 
of this in table 3, in which we compare splittings with and without correlation. To see 
why inclusion of correlation always reduces the splitting, we note the sign and magnitude 
of matrix elements, within ALDA. Even though both (/ 3 ^ LDA ) and (G A C LDA ) are negative, 
because they are dominated by their exchange contributions, we find 

(/ C ALDA ) < -(G ALDA ) < (41) 

because in Eq. (D) antiparallel correlation dominates over parallel correlation. Thus the 
ALDA correlation contribution to the splitting is always negative in SPA. Note that the SIC 



11 



treatment of the splitting is only marginally better than in ALDA because, within SPA, the 
SIC splitting is identical to that of ALDA. 

To analyze the separate levels, we need the magnitude of the SIC corrections: 

(/ ALDA } < (/ o rb) ;= ( H LDA K.,0] } < Q ( (42) 

orika 

but the numerical values of both matrix elements differ by less than 8%. Moreover, 

(G ALDA )>|(/ o rb)|>0 _ (43) 

Looking at the singlet excitation energies of table |I] we see that in ALDA, the s-levels are 
too high (up to 10 mH), whereas the p-levels are too low (by up to 0.4 mH). In X-only 
TDOEP, the s-levels drop (by up to 3 mH), approaching the exact values, but the rise of 
the p-states (by up to 8mH) is too high. Incorporating explicit correlation terms by using 
the TDOEP-SIC kernel, the singlet lines correctly drop further (in comparison to the X-only 
results by up to 1 mH) since 2(/ c ALDA ) - (/ c orb ) « (/ ALDA ) in Eq. (0) is always a negative 
contribution. But still, the p-states are too high. Regarding the triplet excitation energies 
of table |^, the ALDA s-states are too high by at most 6mH, but the p-states are almost 
identical to the exact values. In X-only TDOEP, the triplet states experience a strong 
downshift from the Kohn-Sham excitation energies up to 25 mH, originating from the term 
— (W) (see Eq. (0)). In TDOEP-SIC, this downshift is partly screened by the positive 
correlation contributions 2(G ALDA ) - (/ c orb ), as can be seen from Eqs. (|40|), (@2j) and (|43|). 
This leads to an excellent agreement with the exact values for the s-states. However, these 
correlation terms are too large for the p-states. Since (G ALDA ) > (/ ALDA ), the rise of the 
triplet is always bigger than the dropping of the singlet. 

4.2 Approximate Kohn-Sham potentials 

Next we explore the effect of approximate exchange-correlation potentials v xc on the cal- 
culated excitation spectrum of the He atom. We do not even report results within LDA 
and generalized gradient approximations (GGA)|5(J, since these potentials only support a 
few virtual states, so that many of the transitions reported here do not even exist in such 
calculations. (This problem is worst in small atoms, is less pronounced in molecules, and 
irrelevant in solids). 

To produce a correct Rydberg series, the XC potential must decay as — 1/r, an exact 
exchange effect. Hence we examine the OEP X-only potential (which, for two-electron 
systems is identical to the Hartree-Fock potential) and the OEP-SIC [|3T|] potential. Both 



potentials show the correct behavior for large distances from the nucleus, and support of 
all the Rydberg states is guaranteed. 

Tables |]and ||show the approximate Kohn-Sham excitation energies and the correspond- 
ing corrected excitation energies calculated from the approximate Kohn-Sham eigenvalues 
and orbitals of the X-only potential; tables || and [7] are their analogs from the OEP-SIC 
calculation. The Kohn-Sham orbital energy differences are almost uniformly shifted to larger 
values compared to the orbital energy differences of the exact Kohn-Sham potential. The 
shift ranges from 13.6 mH for the lowest excitation energy to 14.2 mH for excitation en- 
ergies Q n with n > 4 for the X-only potential. The latter shift is exactly the difference 
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Table 4: Singlet excitation energies of neutral helium, calculated from the X-only potential 



and by using approximate XC kernels (in atomic units). 









ALDA (xc) 


TDOEP (X-only) 




Transition 




SPA 


full a 


SPA 


full a 


exact b 


Is -► 


2s 


0.7596 


0.7852 


0.7812 


0.7822 


0.7794 


0.7578 


Is -► 


3s 


0.8533 


0.8598 


0.8601 


0.8588 


0.8591 


0.8425 


Is -► 


4s 


0.8830 


0.8856 


0.8860 


0.8851 


0.8855 


0.8701 


Is -► 


5s 


0.8961 


0.8973 


0.8977 


0.8971 


0.8974 


0.8825 


Is -► 


6s 


0.9030 


0.9037 


0.9040 


0.9036 


0.9038 


0.8892 


Is -» 


2p 


0.7905 


0.7900 


0.7900 


0.7986 


0.7981 


0.7799 


Is -» 


3p 


0.8616 


0.8623 


0.8623 


0.8640 


0.8641 


0.8486 


Is -» 


Ap 


0.8864 


0.8867 


0.8867 


0.8874 


0.8875 


0.8727 


Is -» 


hp 


0.8978 


0.8980 


0.8980 


0.8983 


0.8984 


0.8838 


Is -» 


Qp 


0.9040 


0.9041 


0.9041 


0.9043 


0.9043 


0.8899 


Mean abs 


. dev. c 


0.0118 


0.0156 


0.0153 


0.0162 


0.0161 




Mean percentage error 


1.37% 


1.85% 


1.81% 


1.93% 


1.90% 





a Using the lowest 34 unoccupied orbitals of s and p symmetry, respectively. 

b Nonrelativistic variational calculation [38]. 

c Mean value of the absolute deviations from the exact values. 



Table 5: Triplet excitation energies of neutral helium, calculated from the X-only potential 



and by using approximate XC kernels (in atomic units). 









ALDA (xc) 


TDOEP (X-only) 




Transition 




SPA 


full a 


SPA 


full a 


exact b 


Is 


-► 2s 


0.7596 


0.7493 


0.7488 


0.7370 


0.7345 


0.7285 


Is 


-► 3s 


0.8533 


0.8507 


0.8508 


0.8478 


0.8484 


0.8350 


Is 


^4s 


0.8830 


0.8820 


0.8821 


0.8809 


0.8812 


0.8672 


Is 


-► 5s 


0.8961 


0.8956 


0.8957 


0.8950 


0.8953 


0.8811 


Is 


— ► 6s 


0.9030 


0.9027 


0.9028 


0.9024 


0.9026 


0.8883 


Is 


^2p 


0.7905 


0.7833 


0.7830 


0.7824 


0.7819 


0.7706 


Is 


^3p 


0.8616 


0.8595 


0.8596 


0.8591 


0.8592 


0.8456 


Is 


— > 4p 


0.8864 


0.8855 


0.8856 


0.8853 


0.8854 


0.8714 


Is 


— > 5p 


0.8978 


0.8973 


0.8974 


0.8972 


0.8973 


0.8832 


Is 


— > 6p 


0.9040 


0.9037 


0.9037 


0.9037 


0.9037 


0.8895 


Mean abs. dev. c 


0.0175 


0.0149 


0.0149 


0.0130 


0.0129 




Mean percentage error 


2.11% 


1.78% 


1.78% 


1.53% 


1.51% 





a Using the lowest 34 unoccupied orbitals of s and p symmetry, respectively. 

b Nonrelativistic variational calculation [38]. 

c Mean value of the absolute deviations from the exact values. 
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Table 6: Singlet excitation energies of neutral helium, calculated from the SIC-LDA potential 
and by using approximate XC kernels (in atomic units). 



ALDA (xc) TDOEP (SIC) 



Transition 




SPA 


full a 


SPA 


full a 


exact b 


Is 


-► 2s 


0.7838 


0.8111 


0.8070 


0.8065 


0.8039 


0.7578 


Is 


-► 3s 


0.8825 


0.8891 


0.8895 


0.8878 


0.8881 


0.8425 


Is 


^4s 


0.9130 


0.9156 


0.9161 


0.9150 


0.9154 


0.8701 


Is 


-► 5s 


0.9263 


0.9276 


0.9280 


0.9273 


0.9276 


0.8825 


Is 


-► 6s 


0.9333 


0.9340 


0.9343 


0.9339 


0.9341 


0.8892 


Is 


-> 2p 


0.8144 


0.8145 


0.8144 


0.8222 


0.8217 


0.7799 


Is 


— > 3p 


0.8906 


0.8915 


0.8915 


0.8929 


0.8930 


0.8486 


Is 


— > 4p 


0.9163 


0.9167 


0.9167 


0.9172 


0.9173 


0.8727 


Is 


— > 5p 


0.9280 


0.9282 


0.9282 


0.9285 


0.9285 


0.8838 


Is 


— > 6p 


0.9343 


0.9344 


0.9344 


0.9346 


0.9346 


0.8899 


Mean abs. dev. c 


0.0406 


0.0446 


0.0443 


0.0449 


0.0447 




Mean percentage error 


4.74% 


5.25% 


5.21% 


5.29% 


5.26% 





a Using the lowest 34 unoccupied orbitals of s and p symmetry, respectively. 

b Nonrelativistic variational calculation [38]. 

°Mean value of the absolute deviations from the exact values. 



Table 7: Triplet excitation energies of neutral helium, calculated from the SIC-LDA potential 
and by using approximate XC kernels (in atomic units). 



ALDA (xc) TDOEP (SIC) 



Transition 




SPA 


full a 


SPA 


full a 


exact b 


Is 


-► 2s 


0.7838 


0.7727 


0.7722 


0.7681 


0.7668 


0.7285 


Is 


-► 3s 


0.8825 


0.8799 


0.8800 


0.8786 


0.8789 


0.8350 


Is 


^4s 


0.9130 


0.9120 


0.9121 


0.9115 


0.9117 


0.8672 


Is 


-► 5s 


0.9263 


0.9258 


0.9259 


0.9256 


0.9257 


0.8811 


Is 


— ► 6s 


0.9333 


0.9331 


0.9331 


0.9329 


0.9330 


0.8883 


Is 


^2p 


0.8144 


0.8062 


0.8058 


0.8140 


0.8139 


0.7706 


Is 


^3p 


0.8906 


0.8885 


0.8885 


0.8899 


0.8899 


0.8456 


Is 


— > 4p 


0.9163 


0.9154 


0.9154 


0.9159 


0.9159 


0.8714 


Is 


— > 5p 


0.9280 


0.9275 


0.9276 


0.9278 


0.9278 


0.8832 


Is 


— > 6p 


0.9343 


0.9340 


0.9340 


0.9342 


0.9342 


0.8895 


Mean abs. dev. c 


0.0462 


0.0435 


0.0434 


0.0438 


0.0437 




Mean percentage error 


5.50% 


5.15% 


5.14% 


5.19% 


5.18% 





a Using the lowest 34 unoccupied orbitals of s and p symmetry, respectively. 

b Nonrelativistic variational calculation [38]. 

c Mean value of the absolute deviations from the exact values. 
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Table 8: Singlet-triplet separations in neutral helium calculated from the X-only potential and the 



SIC-potential and by using various approximate XC kernels. Calculated from Eq. (23), using the 
lowest 34 unoccupied orbitals of s and p symmetry All values are in mHartrees. 
x-only SIC 



State 


J xc 


J X 


J xc 


j-TDUKP-SIC 

J xc 


exact 


2S 


32.5 


44.9 


34.9 


37.1 


29.3 


3S 


9.3 


10.7 


9.5 


9.3 


7.4 


AS 


4.0 


4.2 


4.0 


3.7 


2.9 


5S 


2.1 


2.1 


2.1 


1.9 


1.4 


6S 


1.2 


1.2 


1.2 


1.1 


0.8 


2P 


7.1 


16.2 


8.6 


7.8 


9.3 


3P 


2.7 


4.9 


3.0 


3.1 


2.9 


4P 


1.2 


2.1 


1.3 


1.4 


1.3 


5P 


0.6 


1.1 


0.6 


0.7 


0.6 


6P 


0.3 


0.6 


0.4 


0.4 


0.4 


dev. c 


1.46 


4.26 


1.98 


2.26 




rel. dev. 


35% 


49% 


37% 


31% 





including correlation contributions in the form of Vosko, Wilk and Nusair [49] 
b Taken from Ref. @ 

c Mean absolute deviation from the exact values 



between the exact Is eigenvalue (e^ act = —0.90372 a.u.) and the more strongly bound Is 
eigenvalue of the X-only potential (ef,T only = —0.91796 a.u.). Similarly, the Kohn-Sham 
eigenvalue differences calculated in OEP-SIC are shifted by up to 44.5 mH, which again 
is equal to the difference between the Is eigenvalues of the exact Kohn-Sham potential 
and the KS potential in OEP-SIC. In OEP-SIC, the correlation potential is attractive at all 
points in space. Hence, including SIC-correlation contributions into the OEP worsens the 
occupied orbital eigenvalue. To summarize, the inclusion of correlation contributions to the 
ground state potential mostly affects only the occupied state; the virtual states are almost 
exact, i.e., they are almost independent of the choice of the correlation potential. The He 
Kohn-Sham orbitals exhibit a Rydberg-like behavior already for relatively low quantum num- 
bers n f pTXH : already the lower virtual states are mostly determined by the large-r behavior 
of the Kohn-Sham potential, which is governed by the exchange contribution. 

As a consequence, the corrections to the Kohn-Sham orbital energy differences, cal- 
culated on the approximate orbitals, are very close to the corrections calculated from the 
exact Kohn-Sham orbitals. This is most apparent from the singlet-triplet splittings given 
in tables ^|and ||: the splittings depend more strongly on the choice of the XC kernel than 
on the choice of the potential. However, for the excitation energies, the differences among 
the various approximations of the exchange correlation kernel are smaller than the differ- 
ences in the Kohn-Sham excitation energies coming from different potentials. This reflects 
the fact that the resulting orbitals are rather insensitive to different approximations of the 
potential. Hence, the corrections themselves, calculated with approximate XC kernels will 
not cancel the shortcomings of an approximate exchange potential. Tables 0-0 show that 
the corrections go in the right direction only for the singlet states, which are always lower 
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than the corresponding Kohn-Sham orbital energy differences. In other approximations, like 
the LDA and in the popular GGAs for instance, this will be even more severe: There the 
highest occupied orbital eigenvalue is in error by about a factor of two, due to spurious 
self-interaction. There may be error cancellations for the lower Kohn-Sham eigenvalue dif- 
ferences, but in general one should not expect to get a reliable (Kohn-Sham) spectrum in 
LDA and GGAs, because the respective potentials have the wrong behavior for large r. In 
addition, this causes the number of (unoccupied) bound KS states to be finite. 

In total, the inaccuracies introduced by approximate ground state Kohn-Sham potentials 
are substantial, but mostly reside in the occupied eigenvalue for He. It is very unlikely 
that these defects will be cured by better approximations of / xc alone, since the terms 
containing / xc only give corrections to the underlying Kohn-Sham eigenvalue spectrum. 
Hence, the quantitative calculation of excitation energies heavily depends on the accuracy 
of the ground-state potential employed. 



5 Results for the Beryllium Atom 

5.1 Exact Kohn-Sham potential 

The beryllium atom serves as a further standard example for first principles treatments: 
besides numerous quantum chemical studies (e.g. Q), a highly accurate ground-state 



exchange-correlation potential, obtained from quantum Monte-Carlo methods |54j], is avail- 
able for this system. With this potential, we calculated accurate Kohn-Sham orbitals and 
orbital energies of the beryllium atom. In each symmetry class (s, p, and d), up to 38 virtual 
states were calculated on a radial grid similar to the one used in section ^. 

In tables [5] and [ID] we report the excitation energies for the 11 lowest excitations of 
singlet and triplet symmetry. As in helium, the orbital energies of the accurate potential lie 
always in between the experimental singlet and triplet energies. However, the experimentally 
measured singlet-triplet separations in beryllium are much larger than in the helium atom 
(cf. the last columns given in tables |3| and [11]) . Accordingly, to achieve agreement with the 
experimental data, appreciable shifts of the Kohn-Sham eigenvalue differences are needed. 

For the singlet excitation spectrum, given in table |9], the TDDFT corrections yield 
significantly improved excitation energies compared to spectrum of the bare Kohn-Sham 
eigenvalue differences, with average errors reduced by a factor of about 3 regardless of 
which kernel is used. The most distinct improvement towards experiment is achieved for 
the singlet 2P excitation, where the Kohn-Sham eigenvalue difference is off by 32% (61 
mHartree) from the experimental value. 

For the remaining singlet excitations, the TDOEP-SIC kernel yields the best improve- 
ment upon the bare Kohn-Sham spectrum. From figure ^ where the errors for each singlet 
excitation energy are plotted, we see two competing effects: the errors increase with pro- 
gressing angular momentum (with the error of the 3d-states being largest), but decrease 
with progressing principal quantum number n. Note that ALDA has the largest errors for 
the d-states, presumably due to its inability to account for orbital nodes. 

For the triplet spectrum given in table the transition to the 2p state is clearly 
problematic, presumably because of its small magnitude. In particular, the TDOEP X-only 
calculation greatly underestimates the downshift away from the KS eigenvalue difference. 
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Table 9: Singlet excitation energies for the Be atom, calculated from the exact XC potential 
by using approximate XC kernels (in atomic units) 









ALDA 


(xc) 


TDOEP (x-only) 


TDOEP (SIC) 




A; 


-> 3 




CD A 


£,,na 
rull 


DrA 


lull 


QT3 A 

brA 


lull 




is - 


-> Zp 


U.loz i 


U.zU (o 


U.looy 


0.2040 


0.1873 


0.2013 


0.1855 


U.lyoy 


2s 


-> 3s 


0.2444 


0.252b 


n nri r 

0.2515 


n OKIA 
U.z5 l 4 


U.Z05o 


U.z5ot) 


U.z54 / 




2s 


^3p 




u.zoyu 


n 971 A 


0.2748 


0.2758 


0.2739 


0.2750 


fl 97/19 


2s - 


-> 3d 


0.2833 


0.2783 


0.2779 


0.2851 


0.2851 


0.2843 


0.2842 


0.293b 


2s 


-> 4s 


0.2959 


0.2983 


0.2984 


0.2994 


0.2995 


0.2993 


0.2994 


0.2973 


2s - 


-> 4p 


0.3046 


0.3045 


0.3049 


0.3063 


0.3067 


0.3061 


0.3065 


0.3063 


2s ■ 


-» 4d 


0.3098 


0.3084 


0.3084 


0.3106 


0.3106 


0.3104 


0.3103 


0.3134 


2s 


-» 5s 


0.3153 


0.3163 


0.3164 


0.3168 


0.3170 


0.3167 


0.3169 


0.3159 


2s - 


— > 5p 


0.3193 


0.3192 


0.3194 


0.3201 


0.3203 


0.3200 


0.3202 


0.3195 


2s 


-> 6s 


0.3247 


0.3252 


0.3253 


0.3254 


0.3256 


0.3254 


0.3256 


0.325 


2s 


—s- 6p 


0.3269 


0.3268 


0.3269 


0.3273 


0.3274 


0.3272 


0.3273 


0.327 


Mean abs. dev. c 


0.0081 


0.0043 


0.0031 


0.0031 


0.0028 


0.0029 


0.0029 




Mean rel. dev. 


3.75% 


1.69% 


1.15% 


1.27% 


1.09% 


1.13% 


1.15% 




abs. 


dev. d 


0.0028 


0.0033 


0.0029 


0.0025 


0.0025 


0.0024 


0.0024 




rel. dev. d 


0.97% 


1.14% 


1.01% 


0.87% 


0.86% 


0.86% 


0.83% 





a Using the lowest 38 unoccupied orbitals of s, p and d symmetry, respectively. 
b Experimental values from Ref. J55|. 

c Mean value of the absolute deviations from experiment, (all states tabulated) 
Same as c , but excluding the 2s — > 2p transition. 



Because of this effect, we also report average errors with this transition excluded. All 
Kohn-Sham orbital excitations experience a downshift in the ALDA and TDOEP X-only 
calculation. In ALDA, this leads to an overall improvement of the spectrum by a more than a 
factor of 2. The downshift in TDOEP X-only results is too strong, and this behavior is partly 
corrected in the TDOEP-SIC. However, due to overcorrections for the higher excitation 
energies, the average reduction in error over the Kohn-Sham excitation spectrum is only a 
factor of 1.2. 

The errors for the triplet excitation energies are plotted in figure ^|. Clearly, the errors 
of both the Kohn-Sham eigenvalue spectrum and the corresponding corrections decrease 
again with progressing quantum number. Together with the errors plotted in figure ^ this 
signals that the Rydberg-like transitions to states with high principal quantum number n 
are already close to the eigenvalue differences of the accurate Kohn-Sham potential. 

The singlet-triplet separations from equation are given in table |TI| for the three 
different approximate XC kernels. Like in helium, the singlet-triplet splittings are overesti- 
mated by about a factor of two for the S and P transitions if the (diagonal) TDOEP x-only 
kernel is used. The splittings of the D levels, however, appear too small by about a factor 
of two. By the inclusion of correlation contributions to the kernels, the splittings of the S 
and P levels are consistently (and usually correctly) reduced. However, for the D states, 
this correction is always too large, and leads to a reversal of the singlet and triplet energies.^ 

3 This effect can also be observed in the Helium atom. The exact values of the singlet-triplet splittings of 
the D states in helium however, are by two orders of magnitude smaller than in beryllium. 
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Table 10: Triplet excitation energies for the Be atom, calculated from the exact XC potential 
by using approximate XC kernels (in atomic units) 









at rvA 

ALDA 


(xc) 


TDOEP (X-only) 


TDOEP (SIC) 




k - 


3 




CD A 

orA 


£,,na 
rull 


DrA 


An 1 1 a 
IUll 


or A. 


lull 


H/Xpt. 


is — 


-* Zp 


U.loz i 


n nnco 
U.Uyoz 


U.UyUz 


0.0679 


0.0000 


0.0916 


0.0795 


n 1 nno 
U.lUUz 


2s - 


-> 3s 


n O /I A A 

0.2444 


0.2390 


0.2387 


n OQ/in 
U.zo4y 


n OQQQ 

U.zooo 


U.Z4ol 


n 9/1 Qn 


0.2373 


2s - 


^3p 




U.ZOOl 


U.ZOOl 


0.2647 


0.2652 


0.2700 


0.2705 


n 9fi7Q 


2s - 


-> 3d 


r\ o o o o 

0.2833 


0.2807 


0.2805 


0.2814 


0.2813 


0.2867 


0.2865 


n noo'7 

0.2827 


2s - 


-> 4s 


0.2959 


0.2943 


0.2943 


0.2932 


0.2934 


0.2953 


0.2953 


0.2939 


2s - 


-> 4p 


0.3046 


0.3031 


0.3032 


0.3031 


0.3034 


0.3048 


0.3049 


0.3005 


2s - 


-» U 


0.3098 


0.3087 


0.3087 


0.3098 


0.3089 


0.3106 


0.3107 


0.3096 


2s - 


-» 5s 


0.3153 


0.3146 


0.3146 


0.3142 


0.3143 


0.3150 


0.3150 


0.3144 


2s - 


-> hp 


0.3193 


0.3186 


0.3187 


0.3187 


0.3188 


0.3194 


0.3194 


0.3193 


2s - 


-> 6s 


0.3247 


0.3243 


0.3243 


0.3241 


0.3242 


0.3245 


0.3245 


0.3242 


2s - 


-» 6p 


0.3269 


0.3265 


0.3265 


0.3265 


0.3266 


0.3269 


0.3269 


0.3268 


Mean abs. dev. c 


0.0045 


0.0012 


0.0020 


0.0040 


0.0102 


0.0026 


0.0037 




Mean rel. dev. 


3.53% 


0.56% 


1.28% 


3.31% 


9.51% 


1.44% 


2.55% 




abs. dev. d 


0.0017 


0.0012 


0.0012 


0.0012 


0.0013 


0.0020 


0.0020 




rel. dev. d 


0.63% 


0.42% 


0.41% 


0.42% 


0.46% 


0.72% 


0.74% 





a Using the lowest 38 unoccupied orbitals of s, p and d symmetry, respectively. 
b Experimental values from Ref. [55]. 



c Mean value of the absolute deviations from experiment 
d Same as c , but excluding the 2s — > 2p transition. 



From the singlet-triplet splitting in Eq. (|40| ) which, in the SPA, hold for any system since 
the diagonal terms of f XCCTa i cancel, this behavior can be traced back to the overestimation 
of correlation contributions in LDA (in small systems). Self-interaction corrections are not 
expected to cure this shortcoming, for the reason that to leading order the self-interaction 
correction terms cancel in the expressions for the splittings, similar to the way shown in 
section f|. Accordingly, the separations in TDOEP-SIC and ALDA are of similar quality, 
which can be seen from columns one and three in table [TT]. The TDOEP X-only results on 
the other hand, although too small, show the correct ordering of singlet and triplet levels. 

With increasing excitation energy, the difference between the results in SPA and the full 
solution is reduced, as was already observed in the case of helium. The drastic change of 
the triplet 2P state in TDOEP X-only seems to be an artifact of the specific approximation 
to the exchange-correlation kernel, since the results in the SPA and the full calculation for 
this particular excitation energy only differ by 10% if the ALDA is used for f xcaa >. 



5.2 Approximate Kohn-Sham potentials 

The results from using different approximate exchange-correlation potentials for the Be 
atom to calculate the Kohn-Sham eigenvalues and orbitals are given in tables [12] to |15|. 

The errors towards the experimental excitation energies are compiled in figures ^]and |3| for 
the singlet and triplet series. Looking first at the spectra of the bare Kohn-Sham eigenvalues 
(represented in figures ^ and ^| by the points connected with thick lines), we notice that the 
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x-only KLI \ 
x-only KLI v 
x-only KLI \ 
I 



KS Eigenvalue diff . 
ALDA 

TDOEP x-only 

TDOEP SIC 
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ALDA 
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TDOEP x-only 

J I 



4p 

Final State 



6p 



Figure 2: Errors of singlet excitation energies from the ground state of Be, calculated from the accu- 
rate, the OEP-SIC and X-only KLI exchange correlation potential and with different approximations 
for the exchange-correlation kernel (see text). The errors are given in mHartrees. To guide the eye, 
the errors of the discrete excitation energies were connected with lines. 



"HOMO-LUMO" gap is almost independent of the approximation of v xc employed. This is 
in sharp contrast with the He atom case. The correlation contributions cancel for the lowest 
excitation energy, and we must classify this as a non-Rydberg state. For the higher states, 
the situation is different: Starting from the excitation to the 3s level, the series of single- 
particle energy differences appear almost uniformly shifted with respect to the series of the 
exact potential, preserving the typical pattern of their deviation from the experimentally 
measured spectrum. The shifts amount to -14 mH for the OEP-SIC potential, and -34 
mH for the X-only KLI potential. As in helium, these shifts are equal to the differences 
in the eigenvalues of the highest occupied Kohn-Sham orbital: For the accurate potential 
^accurate _ _q.3426 a.u. []56|]), the highest occupied orbitals are more strongly bound than 
in OEP-SIC (e£ Ep - SIC = -0.3285 a.u.) and in X-only KLI (e£r onlyKLI = -0.3089 a.u.). 
Thus, among the virtual states, only the 2p orbital is appreciably influenced by the details 
of the ground state potential. For the higher lying states, the long-range behavior of the 
Kohn-Sham potential dominates. Its —1/r behavior is correctly reproduced both in X-only 
KLI as well as in OEP-SIC. For larger systems, more low-lying excitations can be accurately 
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Figure 3: Errors of triplet excitation energies from the ground state of Be, calculated from the 
accurate, OEP-SIC and X-only KLI exchange correlation potential and with different approximations 
for the exchange-correlation kernel (see text). The errors are given in mHartrees. To guide the eye, 
the errors of the discrete excitation energies were connected with lines. 



approximated, but eventually, for any finite system, the Rydberg excitations will show errors 
due to errors in the ionization potential. Casida et al.[]23[] have studied which excitations 
can be well-approximated with present functional approximations to the potential. 

Regarding the corrections for the singlet excitation energies calculated from Eq. (|23|), 
the first excited state (2p) experiences the largest correction, irrespective of the exchange- 
correlation potential employed. Moreover, the results using different approximate exchange- 
correlation kernels agree within 10 mH. For the remaining singlet excitation energies, the 
calculated corrections using the approximate Kohn Sham orbitals are almost identical to the 
corrections which are obtained from using the accurate Kohn-Sham orbitals. Hence, in figure 
H the errors for the excitations to 3s through 6p show the same pattern of deviations, only 
shifted by the error in the respective eigenvalue of the 2s orbital. On average, the resulting 
singlet excitation energies are closest to experiment, if the approximate exchange-correlation 
potential v xc is combined with the corresponding approximation of the exchange-correlation 
kernel f xca(T >. 

From tables and |15| as well as from figure Bl the behavior of the triplet spectra is 
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Table 11: Singlet-triplet separations in beryllium calculated from the exact XC potential 
by using various approximate XC kernels. Calculated from Eq. (^), using the lowest 38 
unoccupied orbitals of s, p and d symmetry. All values are in mHartrees. 
ALDA TDOEP 



otate 




A-only 




-EjXpt. 


2P 


98.7 


187.3 


106.0 


93.8 


35 


12.8 


21.5 


11.8 


11.8 


3P 


6.2 


10.5 


4.6 


6.4 


3D 


-2.6 


3.8 


-2.3 


10.8 


AS 


4.1 


6.1 


4.1 


3.4 


4P 


1.7 


3.4 


1.6 


5.8 


4D 


-0.2 


1.7 


-0.3 


3.8 


5S 


1.8 


2.7 


1.9 


2.4 


5P 


0.7 


1.5 


0.8 


0.2 


6S 


1.0 


1.4 


1.0 


0.6 


6P 


0.3 


0.8 


0.4 




dev. c 


3.0 


12.4 


3.8 




rel. dev. 


68% 


128% 


75% 





including correlation contributions in the form of Vosko, Wilk and Nusair |49|| . 
b Taken from Ref. gjj 

c Mean absolute deviation from experiment, excluding the 6P state. 



similar, but less unequivocal for the triplet 2p state. For this particular state, the correc- 
tions spread on the order of lOOmH, prevalently due to the significant overcorrection of the 
X-only TDOEP kernels. However, the resulting triplet 2p excitation energy almost exclu- 
sively depends on the approximation to the exchange correlation kernel rather than on the 
exchange-correlation potential employed. On the average, apart from the higher excitations 
in OEP-SIC (c.f. table [14]), the best triplet spectra are obtained if the ALDA is used for 
the exchange-correlation kernel, but this appears to be a fortuitous cancellation of errors. 
The approximate Kohn-Sham excitation energies, except for the 2p state, are already incor- 
rectly lower than the experimental triplet levels. Any further lowering, although correct for 
the eigenvalue-differences of the exact Kohn-Sham potential, actually worsens the triplet 
spectra which are calculated on the basis of an approximate exchange-correlation potential. 
Since the shifts are reduced by correlation contributions in the kernels, the over-corrections 
become less severe for the ALDA and TDOEP-SIC kernels. Another apparent error cancel- 
lation is that when calculating the lowest excitation energy (2s — > 2p) from approximate 
exchange-correlation potentials, the SPA-results are always closer to experiment than the 
full results. This might be related to the fact that for TDOEP X-only, SPA yields the exact 
first-order shift in energy levels in Gorling-Levy perturbation theory|57|], while the "full" 
calculation does not[|]|. In cases where there are large differences between SPA and full 
results, the SPA might be more reliable for these reasons. 

The fact that the corrections to the Kohn-Sham eigenvalue differences only weakly 
depend on the approximation of the exchange-correlation potential v xc , is also reflected 
in table [16[ where the singlet-triplet separations in Be, calculated using the X-only KLI 
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Table 12: Singlet excitation energies for the Be atom, calculated from the X-only KLI 
potential by using approximate XC kernels (in atomic units) 









ALDA (xc) 


TDOEP (X-only) 




k - 


-> 3 




QT3 A 
OrTY 


run 


b-rA 


■Pin 1 1 a 

lull 


h/xpt. 


is — 


-* Zp 


U.lzy / 


0.1990 


0.1795 


0.1958 


0.1791 




2s - 


-> 3s 


U.zloz 


U.ZZ40 


n 99Q9 

U.ZZoZ 


n 99QQ 
U.ZZOO 


U.ZZD 1 


(J.z4yl 


2s - 


^3p 


U.Z4U0 


0.2415 


0.2449 


0.2465 


0.2479 


H 97/19 


2s - 


-> 3d 


0.2527 


n 948D 
U.Z4ioU 




n 9^41 

U.ZO^il 


n 9^ An 


O.z93o 


2s - 


-> 4s 


D 9638 


0.2663 


0.2664 


0.2674 


0.2675 


9973 


2s - 


-> 4p 


0.2725 


0.2727 


0.2735 


0.2745 


0.2751 


0.3063 


2s - 


-» Ad 


0.2773 


0.2758 


0.2759 


0.2780 


0.2780 


0.3134 


2s - 


-» 5s 


0.2822 


0.2833 


0.2834 


0.2838 


0.2840 


0.3159 


2s - 


-> 5p 


0.2863 


0.2864 


0.2867 


0.2872 


0.2876 


0.3195 


2s - 


-> 6s 


0.2913 


0.2918 


0.2919 


0.2921 


0.2923 


0.325 


2s - 


-» 6p 


0.2935 


0.2936 


0.2937 


0.2940 


0.2942 


0.327 


Mean abs. dev. c 


0.0372 


0.0311 


0.0317 


0.0288 


0.0299 




Mean rel. dev. d 


13.5% 


10.4% 


10.7% 


9.5% 


10.1% 




abs. dev. d 


0.0345 


0.0337 


0.0334 


0.0315 


0.0314 




rel. dev. d 


11.49% 


11.18% 


11.07% 


10.39% 


10.38% 





a Using the lowest 38 unoccupied orbitals of s, p and d symmetry, respectively. 

b Experimental values from Ref. J55|. 

c Mean value of the absolute deviations from experiment 

d Same as c , but excluding the 2s — > 2p transition. 



and OEP-SIC potentials are given. The numerical values are close to the results for the 
accurate Be exchange correlation potential in table 0. Again, the obtained splittings are 
more sensitive to the approximation of f xca(T > than to the approximation of the potential 

v xc - 

6 Summary and Conclusion 

In this work we aimed at an assessment of the influence of the three different types of 
approximations (i.e. (i) the XC potential v xc , (ii) the XC kernel / xc and (iii) truncation 
of the space of virtual excitations) inherent in the calculation of excitation energies from 
TDDFT. We calculated the discrete optical spectra of helium and beryllium, two of the 
spectroscopically best known elements, using the exact exchange-correlation potential, the 
KLI-X-only potential and the the KLI-SIC potential for v xc (all three potentials are falling off 
like — 1/r as r — ► oo). These were combined with three approximations for the XC kernel: 
The adiabatic LDA (ALDA), the TDOEP X-only kernel and the TDOEP-SIC kernel. The 
results are given both in the single-pole approximation (SPA) and for a "full" calculation, 
where as many virtual states as possible (typically about 30 to 40) entered the calculation. 
The analysis of these combinations reveals the following trends: First of all, the choice of v xc 
on the calculated spectrum has the largest effect on the calculated spectra. The inaccuracies 
introduced by approximate ground state Kohn-Sham potentials (even those including exact 
exchange) can be quite substantial. This is especially true for the higher excited states, 
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Table 13: Triplet excitation energies for the Be atom, calculated from the X-only KLI 
potential by using approximate XC kernels (in atomic units) 









at rvA 

ALDA 


(xc) 


TDOEP (X-only) 




k - 


-> 3 




CD A 

or A 


lull 


or A 


lull 


Hxpt. 


is — 


-* Zp 


u.izy i 


u.uyoU 


U.UyU ( 


0.0692 


0.0158 


U.lUUz 


2s - 


-> 3s 


U.zlOZ 


r\ o-| i r) 

U.zllz 


n 01 no 
0.2108 


U.zUoy 


U.ZU5 1 


n OQ7Q 

0.23/3 


2s - 


^3p 


U.Z41UD 


U.ZoDZ 


U.ZoDo 


0.2353 


0.2361 


n 9fi7Q 


2s - 


-> 3d 


0.2527 


0.2506 


0.2505 


0.2512 


0.2511 


n noo'7 

0.2827 


2s - 


-> 4s 


0.2638 


0.2622 


0.2622 


0.2611 


0.2613 


0.2939 


2s - 


-> 4p 


0.2725 


0.2710 


0.2711 


0.2709 


0.2712 


0.3005 


2s - 


-» U 


0.2773 


0.2763 


0.2763 


0.2765 


0.2765 


0.3096 


2s - 


-» 5s 


0.2822 


0.2815 


0.2816 


0.2811 


0.2812 


0.3144 


2s - 


-> hp 


0.2863 


0.2856 


0.2857 


0.2856 


0.2857 


0.3193 


2s - 


-> 6s 


0.2913 


0.2909 


0.2909 


0.2907 


0.2908 


0.3242 


2s - 


-» 6p 


0.2935 


0.2931 


0.2932 


0.2931 


0.2932 


0.3268 


Mean abs. dev. c 


0.0300 


0.0291 


0.0298 


0.0323 


0.0371 




Mean rel. dev. 


11.8% 


9.9% 


10.6% 


12.8% 


17.6% 




abs. dev. d 


0.0300 


0.0318 


0.0318 


0.0324 


0.0324 




rel. dev. d 


10.1% 


10.7% 


10.7% 


11.0% 


11.0% 





a Using the lowest 38 unoccupied orbitals of s, p and d symmetry, respectively. 

b Experimental values from Ref. f55| . 

c Mean value of the absolute deviations from experiment 

d Same as c , but excluding the 2s — > 2p transition. 



which appear almost uniformly shifted from the true excitation energies. We observe that 
this shift is closely related to the absolute value of the highest occupied eigenvalue, which, 
in exact DFT, is equal to the first ionization potential of the system at hand. 

For the lower excitation energies, an error cancellation occurs, making these excitations 
less sensitive to the choice of the exchange-correlation potential. This error cancellation 
however, ceases to work the more the excited states behave like Rydberg states. For Helium, 
this is already the case for the first excited state. Hence, in improving the calculation of 
excitation energies from TDDFT requires an improved exchange-correlation potential in the 
first place. The most important requirement for such a potential would be that its highest 
occupied eigenvalue reproduces the experimental ionization potential as closely as possible. 
Empirically, one could introduce a "scissors-operator" similar to the one introduced by Levine 
and Allan [|58|] , which shifts the Rydberg states by a constant being equal to the difference 
between the highest occupied eigenvalue and the negative of the experimental ionization 
potential. But such a procedure would not produce a first-principles calculation. In our 
opinion however, the construction of approximate exchange-correlation potentials based on 
orbital functionals would be the method of choice for the future. 

The effect of the choice of the exchange-correlation kernel on the calculated spectra, 
in turn, is much less pronounced. However, its relative importance increases whenever the 
"first-order" effects, originating from v xc cancel. This is the case for the values of the 
singlet-triplet splittings and the lower excitation energies of Be. For these "second-order 
effects", the correlation contributions contained in f xc are important. We observe that the 
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Table 14: Singlet excitation energies for the Be atom, calculated from the SIC-KLI potential 
by using approximate XC kernels (in atomic units) 
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n 9£ 1 a 




97/19 
U.Z ( 4Z 


2s - 


-> 3d 


0.2684 


0.2632 


0.2628 


9694 
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0.2936 


2s - 


-> 4s 


D 9890 


9844 


9846 


0.2855 


0.2856 


9973 


2s - 


-> 4p 


0.2909 


0.2910 


0.2916 


0.2926 


0.2931 


0.3063 


2s - 


-» U 


0.2956 


0.2943 


0.2943 


0.2961 


0.2962 


0.3134 


2s - 


-» 5s 


0.3013 


0.3023 


0.3025 


0.3028 


0.3030 


0.3159 


2s - 


-> hp 


0.3054 


0.3054 


0.3056 


0.3062 


0.3064 


0.3195 


2s - 


-> 6s 


0.3106 


0.3112 


0.3113 


0.3114 


0.3116 


0.325 


2s - 


-» 6p 


0.3129 


0.3129 


0.3130 


0.3133 


0.3134 


0.327 


Mean abs. dev. c 


0.0210 


0.0155 


0.0153 


0.0130 


0.0138 




Mean rel. dev. d 


8.07% 


5.29% 


5.25% 


4.32% 


4.78% 




abs. dev. d 


0.0168 


0.0161 


0.0158 


0.0140 


0.0139 




rel. dev. d 


5.65% 


5.35% 


5.26% 


4.60% 


4.60% 





a Using the lowest 38 unoccupied orbitals of s, p and d symmetry, respectively. 

b Experimental values from Ref. J55|. 

c Mean value of the absolute deviations from experiment 

d Same as c , but excluding the 2s — > 2p transition. 

ALDA for the XC kernels already leads to quite reasonable results which are only marginally 
improved by using the more complicated TDOEP-SIC kernel. Besides the missing frequency 
dependence, correlation contributions are hard to model on top of an exact exchange treat- 
ment and one might speculate that the ALDA takes advantage from a fortuitous error 
cancellation between exchange and correlation effects. Again, we expect only orbital func- 
tional to manage a marked improvement over the ALDA, which, up to now, can has been 
the workhorse of TDDFT. 

Finally, the inevitable truncation of the space of virtual excitations is appreciable only for 
the lowest lying states. In most cases, the results of the single-pole approximation (SPA), 
which, in the nondegenrate case, merely requires a pair of "initial" and "final" KS states 
are close to the results obtained from using more configurations. 

We thank C. Umrigar for providing us with the data of the exact exchange-correlation 
potentials for helium and beryllium. This work was supported in part by the Deutsche 
Forschungsgemeinschaft and by NATO. K.B. acknowledges support of the Petroleum Re- 
search Fund and of NSF grant CHE-9875091. 
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Table 15: Triplet excitation energies for the Be atom, calculated from the SIC-KLI potential 
by using approximate XC kernels (in atomic units) 
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c Mean value of the absolute deviations from experiment 
d Same as c , but excluding the 2s — > 2p transition. 
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